/*
[] brackets indicate internal census variable name that can't be disclosed. 

Randomly select piks and store the seeds from each run through the loop. 
1) Create an empty base table on which to append the seeds (only do this once before running). 
2) Pull in the seed as a macro variable using the loop counter from the bash script.
3) Randomly sample using proc surveyselect, print seed output, and append seed to base table. 
4) Merge back to original data for regressions. 
*/

libname hannah '/projects/users/########/Snapshot2022/IntermediateData';
libname hantmp '/projects/users/########/Snapshot2022/IntermediateData/TempData';
libname ans '/projects/users/########/Snapshot2022/AnalysisData';

/*
* initialize a base empty seed table for appending later seeds;
* ONLY RUN THIS ONCE - DON'T OVERWRITE ON EACH RUN THROUGH BASH LOOP;
proc sql;
  create table hantmp.pikloopseeds
    (
    Label1 char(21) format=$21. informat=$21. label='Label1',
    cValue1 char(13) format=$13. informat=$13. label='cValue1',
    nValue1 num format=D12.3 label='nValue'
    );
quit;
*/

/* 
* get counties merged to individual piks (could do hex IDs);
* ONLY NEED TO DO THIS ONCE;
proc sql;
  create table hantmp.nohigh_pikscnt as select
  a.*, b.*
  from hantmp.rand_pikscnt as a 
    inner join hantmp.pikspine as b on a.pik=b.pik
  order by a.pik;
quit;

data hantmp.nohigh_pikscnt;
  set hantmp.nohigh_pikscnt (keep = pik [county id] nohigh);
  if nohigh=0 then delete;
  drop nohigh;
run;

* drop duplicates - get unique piks and their county IDs;
proc sort data = hantmp.nohigh_pikscnt nodupkey;
  by pik;
run;

* data must be sorted by strata;
proc sort data = hantmp.nohigh_pikscnt;
  by [county id];
run;
*/

/* 
* Already did - ONLY NEEDED TO DO THIS ONCE IF I STORE THESE COUNTY NSIZE;
* get frequencies in each strata;
proc freq data = hantmp.nohigh_pikscnt;
  tables [county id] / out=nohigh_2freq;
run;

* determine strata sizes. want 1000000 obs. need to define a dataset with that number;
*_nsize_ is the sample size from each strata, positive integer;
* If not a positive integer, deleted and error dataset created;
data nohigh_2freq2 nohigh_2error;
  set nohigh_2freq;
  sampnum=(percent * 1000000)/100;
  _nsize_ = round(sampnum,1);
  sampnum = round(sampnum,0.01);
  if _nsize_ = 0 then output nohigh_2error;
  if _nsize_ = 0 then delete;
  if count >= ##### then _nsize_ = #####;
  if count < ##### then _nsize_ = count;
output nohigh_2freq2;
run;
* keep strata variable and sample size in each stratum;
data nohigh_2freq3;
  set nohigh_2freq2;
  keep [county id] _nsize_;
run;
data nohigh_2err3;
  set nohigh_2error;
  keep [county id] _nsize_;
run;
proc append base = nohigh_2freq3 data = nohigh_2err3;
run;
proc sort data = nohigh_2freq3;
  by [county id];
run;
*/

* Pull in seed from bash loop counter and save as a macro variable;
%let seedi=%sysget(VAR1);

* randomly select one million piks and append the seed to the base table;
proc surveyselect data = hantmp.nohigh_pikscnt
  out = hantmp.nohigh_pikloop_lpdid
  sampsize = hantmp.nohigh_2freq3
  seed = &seedi;
  strata [county id];
run;

* now need to merge back to full data to make a file Stata can analyze;
proc sql;
  create table ans.nohigh_pikloop_lpdid as select
  a.*, b.*
  from ans.nohigh_analysis as a 
    inner join hantmp.nohigh_pikloop_lpdid as b on a.pik=b.pik
  order by a.pik, a.year;
quit;


/* 
*For specific groups (e.g., nohigh);
ods output Summary = Summary;
proc surveyselect data = hantmp.nohigh_piks
  out = hantmp.nohigh_pikloop
  method = srs
  sampsize = 1000000
  seed = config outseed;
run;

proc append base = hantmp.pikloopseeds data = summary;
run; quit;
*/

